Physical Pictures of Transport in Heterogeneous Media: Advection-Dispersion, 
Random Walk and Fractional Derivative Formulations 
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The basic conceptual picture and theoretical basis for development of transport equations in 
porous media are examined. The general form of the governing equations is derived for conservative 
chemical transport in heterogeneous geological formations, for single realizations and for ensemble 
averages of the domain. The application of these transport equations is focused on accounting for 
the appearance of non-Fickian (anomalous) transport behavior. The general ensemble-averaged 
transport equation is shown to be equivalent to a continuous time random walk (CTRW) and re- 
duces to the conventional forms of the advection-dispersion equation (ADE) under highly restrictive 
conditions. Fractional derivative formulations of the transport equations, both temporal and spa- 
tial, emerge as special cases of the CTRW. In particular, the use in this context of Levy flights is 
critically examined. In order to determine chemical transport in field-scale situations, the CTRW 
approach is generalized to non-stationary systems. We outline a practical numerical scheme, sim- 
ilar to those used with extended geological models, to account for the often important effects of 
unresolved heterogeneities. 



1. Introduction 

Quantification of chemical transport mediated by flow 
fields in strongly heterogeneous geological environments 
has received an inordinate amount of attention over the 
last three decades, and a vast literature dealing with the 
subject has developed (see, e.g., the recent reviews in Da- 
gan and Neuman [1997]). Existing modeling approaches 
are generally based on various deterministic and stochas- 
tic forms of the advection-dispersion equation (ADE); 
the former include conditioning the domain of interest 
by known heterogeneity structures, while the latter in- 
clude Monte Carlo, perturbation and spectral analyses. 
A major feature of transport, particularly in more hetero- 
geneous domains, is the appearance of "scale-dependent 
dispersion" [e.g., Gelhar et at, 1992]. Contrary to the 
fundamental assumptions underlying use of the classical 
ADE (which assumes a constant flow field and dispersion 
coefficients) , the very nature of the dispersive transport 
seems to change as a function of time or distance traveled 
by the contaminant. Such scale-dependent behavior, also 
sometimes referred to as "pre-asymptotic" , "anomalous" 
or "non-Gaussian" , is what we shall refer to as "non- 
Fickian" transport. 

Efforts to quantify non-Fickian transport have focused 
on more general stochastic ADE's with, e.g., spatially 
varying velocity fields. Stochastic analyses have provided 
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substantial insight into the dispersion process. They have 
been shown, through application to well-documented 
field experiments, to provide predictions of the tempo- 
ral variation of the first and second order moments of 
tracer plumes in geological formations characterized by 
relatively small degrees of heterogeneity (e.g., the Cape 
Cod site [Garabedian et al, 1991]). Other variations 
based on the classical ADE have also received attention; 
these include "patch" solutions which include an empir- 
ical time- or space-dependent dispersivity, and mobile- 
immobile and multirate diffusion type models [e.g., Hag- 
gerty and Gorelick, 1995; Harvey and Gorelick, 2000]. 
However, the vast majority of these models assume, ei- 
ther explicitly or implicitly, an underlying Fickian trans- 
port behavior at some scale [e.g., Sposito et al., 1986; 
Rubin, 1997]. Also, many of these approaches are based 
on perturbation theory, and they are therefore limited to 
porous media in which the variance of the log hydraulic 
conductivity is small. 

Other non-local formulations that do not invoke a Fick- 
ian transport assumption have been hypothesized and/or 
developed from various mathematical formalisms [e.g., 
Zhang, 1992; Glimm et al, 1993; Neuman, 1993; Deng 
et al., 1993; Cushman et al., 1994; Dagan, 1997]. These 
formalisms, in general, are founded on a fundamental sep- 
aration between advective and dispersive mechanisms; 
they yield solutions (for the concentration) that result 
in definition of a dispersion tensor that is usually formu- 
lated in Fourier-Laplace space, whose inversion is difficult 
to treat and/or apply. 

Practical application of these models, to quantify the 
full evolution of a migrating contaminant plume, has not 
yet been achieved. In fact, the overwhelming emphasis of 
these various studies has been limited to moment char- 
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acterizations of tracer plume migration, and/or to deter- 
mination of the "macrodispcrsion" parameter. The com- 
plete solutions are not analytically tractable, and their 
practical utility remains largely undemonstrated. 

The difficulty in capturing the complexities of tracer 
plume migration patterns suggests that local, small-scale 
heterogeneities cannot be neglected. Evidently, these un- 
resolvable heterogeneities contribute significantly to the 
occurrence of non-Fickian transport. The apparent ex- 
istence of hydraulic conductivity fields with coherence 
lengths that vary over many scales suggests that tempo- 
ral, as well as spatial issues must be considered in any 
mathematical formulation. Coupled to this problem is 
the lack of clarity of how best to use field observations 
to reduce the inevitable uncertainties of the model. Fre- 
quently, the latter issue involves the interplay between 
ensemble averaging (probabilistic approaches) and spa- 
tial scales of resolution of non-stationary geological fea- 
tures. 

In this paper, we re-evaluate the basic conceptual pic- 
ture of tracer migration in heterogeneous media. We de- 
rive the general form of the governing equations for con- 
servative chemical transport in heterogeneous geological 
formations, for single realizations and for ensemble aver- 
ages of the domain. We emphasize quantification of non- 
Fickian transport behavior, and show that a general form 
of the ensemble-averaged transport equation is a contin- 
uous time random walk (CTRW). In this framework, we 
show that non-Fickian transport results from the inap- 
plicability of the central limit theorem to capture the 
distribution of particle transitions (detailed in the next 
section). Fractional derivative formulations of the trans- 
port equations, both temporal and spatial, are seen to 
emerge from another set of conditions, and are therefore 
special cases of the CTRW. We then focus on quantify- 
ing transport in non-stationary media, and discuss how 
best to deal with the coupled problem of integrating en- 
semble averaging with information on non-stationarity at 
various scales of resolution. 



2. Governing Transport Equations for 
Heterogeneous Media 

2.1. Physical Framework of the Transport Equations 

Contaminants disperse as they migrate within the flow 
field of the geological maze we call an aquifer. At the out- 
set one must choose an underlying physical model of this 
process. Two possible models include Taylor dispersion 
and multiple transitions. Taylor dispersion is based on 
molecular diffusion of particles in a flowing fluid (e.g., 
in a pipe) and is governed by an ADE, to be discussed 
below. An identical formulation can be obtained by con- 
sidering particle movement in a random network and ap- 
plying the central limit theorem. The extensive use of the 
ADE in the hydrology literature is based essentially on 
the generic concept of Taylor dispersion and works well 



for relatively homogeneous systems. The particles are as- 
sumed to be transported by the average flowing fluid in 
the medium while the "diffusion" is the dispersion due 
to local medium irregularities. Larger scale effects (e.g., 
permeability changes) are treated as perturbations of this 
model in conventional stochastic treatments. 

The prime interest in this work is in highly heteroge- 
neous systems; in these systems contaminant motion can 
be envisioned as a migrating cloud of particles, each of 
which executes a series of steps or transitions between 
changes in velocity v. The spatial extent of these tran- 
sitions depends on the criterion used to define changes 
in v. The classical approach is to consider the system 
divided into representative elementary volumes (REVs) 
and determine an average v and dispersion D in each 
REV. In our approach we dispense with the REV idea, 
because averages can be unreliable in a system of very 
wide fluctuations about the mean value. The change of 
concentration AC at each position in a time increment 
At is At x (the net particle flux). The effective volume 
contributing the net particle flux in At can vary con- 
siderably at different positions in the system. Thus the 
length scale over which AC varies slowly in space can 
change considerably over the system. If one fixes a sam- 
pling volume at each position, it is important to retain 
the full distribution (not an average) of the transition 
times (determined with a physical model) of flux con- 
tributing to AC If this distribution is retained, then 
in our approach one can still use the limit of a spatial 
continuum (as shown below). 

The distribution of transition times, ip(t), can be deter- 
mined in principle from an analysis of the streamtubes 
of the flow field and contains the subtle features that 
can produce non-Fickian behavior. The physical features 
necessary for non-Fickian transport are the existence of 
a wide range of transition times (causing large differ- 
ences in the flow paths of migrating particles) and suffi- 
cient encounter with statistically rare, but rate-limiting 
slow transitions (e.g., low velocity regions) [Berkowitz 
and Scher, 1995]. These general ideas will be developed 
schematically in the next sections. 



2.2. Single Realization Transport Equation 

For our point of departure we need a transport equa- 
tion framework that can enumerate all these possible 
paths and encompass the motion from continuous to dis- 
crete over a range of spatial and temporal scales, for 
any given realization of the domain. An excellent candi- 
date is the "Master Equation" [Oppenheim et al., 1977; 
Shlesinger, 1996] 



= - ]T w(s', s)C(s, t) + ^ u;(s, s')C(s', t) 

s' s' 

(1) 
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for C(s,t), the particle concentration at point s and time 
t, where w(s, s') is the transition rate from s' to s (the di- 
mension of T, s w is reciprocal time). The transition rates 
describe the effects of the velocity field on the particle 
motion; the determination of w(s, s') involves a detailed 
knowledge of the system. We assume the average effective 
range of w(s, s') is a finite distance. The Master Equa- 
tion has been applied in the context of electron hopping 



in random systems [e.g., Klafter and Silbey, 1980a], and is 
discussed widely in the physics and chemistry literature. 

The transport equation in (|l|) does not separate the 
effects of the varying velocity field into an advective and 
dispersive part of the motion; this separation is an ap- 
proximation based on the assumption of relatively homo- 
geneous regions in which C(s, t) will be slowly varying 
over a finite length scale (the range of transition rates), 



C(s', t) w C(s, t) + (s' - s) • VC(s, t) + \ (s' - s) (s' - s) : VVC(s, t) (2) 

(with the dyadic symbol : denoting a tensor product). Substituting (0) into (R) leads to a continuum description (i.e., 
local diffusion in a pressure field 7r(s)) and a partial differential equation (pde), for a single realization of the domain: 

= ^Ws, S ')-4', S ))C(M)+^^ lS ')(s'-s).VC(M) 

s' s' 

+ ^w(s,s')i(s' - s)(s' - s) : VVC(s,i). (3) 



r 



We note that (||) is close to the form of an ADE with 
the exception of the term proportional to C(s,t). This 
term is present due to the asymmetry of the transition 
rates (due to the bias of the pressure field) and/or the 
non-stationary medium (due to the explicit position de- 
pendence of the rates - cf . (^)). It makes a contribution 
to the final form of the pde for diffusion in a force field. 
If the system is stationary this term vanishes (as we show 
below) and thus reduces to the form of an ADE. One can 
already observe in (j|) generalized velocity and dispersion 
coefficients (in terms of w(s, s')); however we have not yet 
separated out the effects of the flow field and determined 
transport coefficients. In order to fully determine the 
final pde and separate the advection and diffusion con- 
tributions, we must specify the w(s,s') in terms of 7r(s), 
the pressure field. 

A general form for a non-stationary medium is 

w(s, s') = W (|s' - s|; s') Q(tt(s') - tt(s)) (4) 

where the asymmetry in the rates is due to 7r(s') — 7r(s), 
the pressure difference at s' and s, and the explicit depen- 
dence of the overall rate W on location (Q is a function of 
the pressure difference only) . We specify the fi-function, 
so that @ is written as 

w(s,s') = W(\s' -s|;s')0(7r(s / ) -tt(s)) 



» F(|s'- S |; S ')[A+i(7r(s')-7r( S ))] (5) 

where in (^|) non-linear terms in the pressure difference 
have been neglected (i.e., terms proportional to (V7r) 2 ) 
and a contribution to the transition rates is retained even 
for vanishing pressure difference. The significance of the 
latter step can be seen by realizing that F(tt(s')—it(s)) is 
a simple advection contribution (with a permeability pro- 
portional to F) and the term FA is proportional to a local 
diffusion contribution to the rates. The A term retains 
the scattering effects of the medium (i.e., the transfers 
between "streamtubes" ) even in the limit of very small 
local pressure differences. It is also closely associated 
with the effect of "local" dispersion. 

We now also assume F(|s' — s|; s') will be slowly varying 
over some finite length scale. We expand in a Taylor 
series to second order in s' — s, 

F(|s'-s|;s') ps F(|s'-s|;s) + (s'-s)- VF 

+§(s'-s)(s'-s) : VVF. (6) 

In (0), the gradient operates on the second argument, s'. 
Combining (^|) and (g), and substituting into the first 
term on the right side of (^), we have 



w(s, s') - w(s', s) « F(\s' - s|; s') [A + ±(tt(s') - w(s))] - F(\s' - s|; s) [A - ±(7r(s') - 7r(s))] 
« F(|s'-s|; S )(4s')-^(s)) 

+ [(s' - s) ■ VF + i(s' - s)(s' - s) : VVF] x [A + i(7r(s') - tt(s))] . (7) 
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Now using a similar expansion for the pressure difference, 
we have 



7r(s')-7r(s) ps (s' - s)-Vtt(s)+±(s' - s)(s' - s) : VVtt(s) 
Substituting (||) into (0) and using 

X>(|b'- S |;s)(s'-s)=0 



(8) 
(9) 



because F is an even function of the vector difference, we 
obtain for the expression in (Q), summed over s', 



D(s) 



Vtt(s) + V • VD(s) 



(10) 



where the dispersion tensor is defined as 



D(s) = ± £ F(|s' - s|; s)(s' - s)(s' - s)A. (11) 

s' 

We insert @ into (§) and use (§-(§), (§), (|) (cf. Ap- 
pendix A) to obtain 



dC{B,t) 

at 



D(b) 
A 



Vtt(s)C(s, i) + V(D(s)C(s, t)) 



(12) 

The form of ([12]) is a continuity equation - the time 
derivative of the concentration is equal to the divergence 
of the total concentration flux, the sum of the diffusive 
concentration flux and the advective concentration flux 
- with an effective permeability of 



k(s) 



D(s) 
A ' 



(13) 



Equation ( Jl 2| ) is a generalization to a non-stationary 
medium of the well-known Smoluchowski equation 
[Chandrasekhar, 1943] which is the basis for describing 
diffusion in a force field. In our case the force field is 
Vtt(s). In the case of electron transfer in a potential 
field the A in ([l3]) can be shown to be kT (where T is 
the temperature and k is Boltzmann's constant) and the 
relation in ( |l3| ) is the Einstein relation between mobil- 
ity and diffusion. We use a convention that a product 
between a tensor T and a vector V is TV yielding a vec- 
tor. In our case, the vector v(s) = — k(s)V7r(s) is the 
velocity field and for an incompressible fluid V-v(s) = 0. 
The only term remaining in (12) proportional to C(s,i) 
is V • VD(s)C(s,i). The final form for the pde for an 
incompressible fluid is 



dC(s,t) 
dt 



-v(s) • VC*(s, t) + V • V(D(s)C(s, t)). (14) 



Equation (Q) is a generalization of the ADE. While 
many simplifications of the ADE are based on (14) with 



D(s) = D (i.e., a constant), the usual ("general") form of 
the ADE includes a s-dependent D in ( |l4| ) but with the 
second term replaced by V ■ (D(s)VC(s, *)). Thus 
differs from this usual form of the ADE by the addition of 
two terms: V • VD(s)C(s,i) and VD(s) • VC. The form 
of (|l4| ) is the same as postulated by Kinzelbach [1986], 
based on the Ito process. 

The difference in the general form of the ADE can be 
traced to starting the derivation with the pressure field 
7r(s) and not with V7r(s), i.e., the expansion (ph is treated 
on the same basis as the other expansions J0) and (|^). 
Hence, starting with the Master equation (|l]j and using 
a general expression for the transfer rates we obtain, for 
a specific heterogeneous medium, in a continuum limit 
(slowly varying C(s, t) and w(s, s')) the generalized equa- 
tion for diffusion in a force field (Smoluchowski) which 
for irrotational flow is a generalized ADE. We assert that 
for a non-stationary medium, i.e., s-dependent v and D, 
(|l4| ) should be the starting point for numerical calcu- 
lations. The main numerical differences between this 
equation and the usual ADE (with D(s)) should arise 
in "boundary" regions of more spatially varying D(s). 
The importance of accounting for D(s) has been demon- 
strated by, e.g., Labolle et al. [1996]. 

We will show that the "standard" ADE emerges as the 
continuum limit of the ensemble averaged Master equa- 
tion (the term proportional to C(s, t) vanishes for sta- 
tionary transition rates) . In general, the continuum limit 
presents difficulties in regions of increased heterogeneity, 
such as tightly interspersed permeability layers. The con- 
centration C(s,t) will not necessarily vary slowly on the 
same length scale throughout the system. The point av- 
erage of v and D can be very sensitive to small changes 
in the local volume used to determine the average. Con- 
versely, if one fixes the volume to a practical pixel size 
(e.g., 10 m 3 ) the use of a local average v and D in each 
volume can be quite limited, i.e., the spreading effects of 
unresolved residual heterogeneities are suppressed [e.g., 
Dagan, 1997]. We will return to this issue in a broader 
context in section 4. It essentially involves the degrees of 
uncertainty and its associated spatial scales. We start, at 
first, with an ensemble average of the entire medium and 
discuss the role of this approach in the broader context. 



2.3. Ensemble Average Transport Equation 

We resume our examination of the Master Equation 
approach, i.e., before assuming any continuum limit. The 
ensemble average of (|l|) can be shown [Klafter and Silbey, 
1980b] to be of the form 
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dP(s,t) 

at 



V / cf>(s' -s,t-t')P(s,t')dt' + J2 [ <p(s-s',t-t')P(s',t')dt' 
a , Jo a , Jo 



(15) 



where P(s,t) is the normalized concentration, and </>(s,i) 
is defined below in (^Oj) . The form of ([l5|) is a "General- 
ized Master Equation" (GME) which, in contrast to (jl]), 
is non-local in time and the transition rates are station- 
ary (i.e., depend only on the difference s - s') and time- 
dependent. This equation describes a semi-Markovian 
process (Markovian in space, but not in time), which ac- 
counts for the time correlations (or "memory" ) in particle 
transitions. 

It is straightforward to show [Kenkre et al., 1973; 
Shlesinger, 1974], using the Laplace transform, that the 
GME is completely equivalent to a continuous time ran- 
dom walk (CTRW) 



„, Jo 



,t~t')R(s',t')dt' (16) 



where i?(s,i) is the probability per time for a walker to 
just arrive at site s at time t, and i/j(s, t) is the probability 
rate for a displacement s with a difference of arrival times 
of t. The initial condition for P(s, t) is S s ,o3(t— + ), which 
can be appended to (fl6|). The correspondence between 
(|l|) and |l|) is 



P(s,i) = [ <H(t-t')R(s,t')dt' 
Jo 



where 



(i) = 1 - / ip(t')dt' 
Jo 

is the probability for a walker to remain on a site, 



^(t)=$^(M) 



and 



(s,u) 



utp(s, u) 
1 - ip(u) 



(17) 



(18) 



(19) 



(20) 



where the Laplace transform (£) of a function f(t) is 
denoted by f(u). 



Equations ©-© are in the form of a convolution 
in space and time and can therefore be solved by use of 
Fourier and Laplace transforms [Scher and Lax, 1973]. 
The general solution is 



P(k,«) 



i - 4>(u) 



1 - A(k,u) 



(21) 



where P(k, u), A(k, u) are the Fourier transforms (F) of 
P(s,u), -0(s,u), respectively. 

The CTRW accounts naturally for the cumulative ef- 
fects of a sequence of transitions. The challenge is to 
map the important aspects of the particle motion in the 
medium onto a ifj{ s jt)- The identification of tp( s it) nes 
at the heart of the CTRW formulation. The CTRW ap- 
proach allows a determination of the evolution of the par- 
ticle distribution (plume), P(s,t), for a general iji(s,t); 
there is no a priori need to consider only the moments 
of P(s^t). As we discuss below, a ip(s,t) with a power 
law (pG) for large time leads to the description of anoma- 
lous transport (e.g., non-Fickian plumes). Once V , ( s ^) 
is defined one needs to calculate A(k, u) and then de- 
termine the propagator P(s,i) by inverting the Fourier 
and Laplace transform of (pl|). The latter can be quite 
challenging. 

As shown previously the separation between advection 
and dispersion occurs in the continuum (diffusion) limit. 
In an ensemble averaged system this limit leads to an 
ADE [Berkowitz and Scher, 2001]. For clarity and con- 
venience, we reproduce the argument here. The first step 
is to make a series expansion of P(s, t) similar to (||); in- 
serting this into (|l^) yields 



8P{s,t) 
dt 



./ Jo 



(s -s',t- t')(s' - s) • VP(s, t') + 4>(s -s',t- t')±(s' - s)(s' - s) : VVP(s, t')]. 



(22) 



We write (22) in a more compact form 



dP(s,t) 



= f tft'[-v v ,(t-t')-VP(s,0 + *v>(*-0: VVP(s,0] 
Jo 



(23) 
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(24) 



(25) 



r 



Note the sum (over s') in ( |22| ) is independent of s in a 
stationary system; hence we shift the summation vari- 
able to obtain (p4|)-(|25|). This particular formulation is 
convenient because, in (|23|), we can define terms that are 
familiar in the context of traditional modeling: the "effec- 
tive velocity" v^, and the "dispersion tensor" <&,/,. Note, 
however, that both of these terms are time-dependent, 
and most significantly, depend fundamentally on ip(s,t). 
This equation has the form of an ADE generalized to 
non-local time responses as a result of the ensemble av- 
erage. 

The next step is a crucial one in distinguishing between 
normal and anomalous transport. If ^>(s, t) has both a 
finite first and second moment in t the transport is nor- 
mal and one can expand ip(s,u) as [Scher and Montroll, 
1975] 



i/j(s,it) = pi(s) - p 2 (s)u + p 3 (s)u 2 + 
and ip(u) = -0(s, u) = 1 — iu + du 2 - 



(26) 



with X) S P 1 ( S ) = 1' the normalization of ^(s,i), and 
J2 S P 2 ( S ) = t an d J2 S P3( S ) = d, the nrs t an d second 
temporal moments of ijj(t), respectively. Note that small 
u corresponds to large time in Laplace space. The func- 
tions Pi(s) are asymmetric due to the bias in the velocity 
field; pi(s) is the probability to make a step of displace- 
ment s. One now inserts (£6|) into ( |20| ) and expands in 
a power series of u. The leading term is independent of 
u, which we retain. The correction to this leading term 
is proportional to u and is small. Substituting this ex- 
pression into the Laplace transform of (p3|)-(|25|), which 
is (|5^)-([35]) (cf. below), and taking the inverse Laplace 
transform of the result, yields the ADE 



dP{s,t) 
dt 



-v ■ VP(s,t) +D : VVP(s,t) 



(27) 



where the effective velocity v is equal to the first spatial 
moment of pi(s),s, the mean displacement for a single 
transition, divided by the mean transition time i, and the 
dispersion tensor D = is the second spatial moment 
divided by t, which can be written as 



^pi(s)s/i = s/t 



(28) 



(29) 



where v = |v| and s = |s|. If we retain the term propor- 
tional to u when inserting ( p6| ) into ( p0| ) , we obtain terms 
with both spatial and temporal derivatives of P(s,t). 



Thus, our underlying physical picture of advective- 
driven dispersion reduces to the familiar ADE when one 
can assume smooth spatial variation of P(s, t) and finite 
first and second temporal moments of i(j(s, t). 



2.4. Non-Fickian Dispersion 

When the ip(s,t) has a power law (algebraic tail) de- 
pendence on time at large t, i.e., 



#s,t)~f 



(30) 



the first and second temporal moments do not exist for 
< (3 < 1, while the second temporal moment does not 
exist for 1 < (3 < 2. The dependence of ip(s,t) in ( |30| ) is 
a manifestation of a wide distribution of event times as 
encountered in highly heterogeneous media. The relation 
between the power law behavior (|3^) and non-Fickian 
(anomalous) transport has been well documented [e.g., 
Scher and Montroll, 1975; Berkowitz and Scher, 2001]. 
We sketch the key points of that relationship: The form 
of ip(s,t) at large time determines the time dependence 
of the mean position £(t) and standard deviation a(t) of 
P(s, t). In the presence of a pressure gradient (or "bias" ), 
and for (p(i|), it can be shown [Scher and Montroll, 1975; 
Shlesinger, 1974] for < (3 < 1 that 







while for 1 < (3 < 2 



e(t) - t 



£{t) - t 



(31) 



(32) 



(33) 



(34) 



Moreover, it can be shown that Fickian-like transport 
arises when (3 > 2 [e.g., Margolin and Berkowitz, 2000] . 

The unusual time dependence of i(t) and a(t) in (pl|)- 
([m]), resulting from the infinite temporal moments of 
ip(s,t) (i.e., the conditions of the central limit theo- 
rem are not fulfilled) , is the hallmark of the non-Fickian 
propagation of P(s,t). This behavior is in sharp con- 
trast to Fickian models where, I(t) ~ t and a(t) ~ t 1 ^ 2 
(as an outcome of the central limit theorem) and the 
position of the peak of the distribution coincides with 
£(i). Note that in Fickian transport, £(t)/a(t) - t 1 / 2 ; 
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an important distinguishing feature of anomalous trans- 
port is that £(t)/a(t) ~ constant for < f3 < 1, and 
l(t)/a(i) ~ iO 3 " 1 )/ 2 for 1 < /3 < 2. The relative shapes 
of the anomalous transport curves, and the rate of ad- 
vance of the peak, vary strongly as a function of p. Thus 
the parameter P effectively quantifies the contaminant 
dispersion; this parameter is discussed in detail by, e.g., 
Margolin and Berkowitz [2000, 2002] and Berkowitz and 
Scher [2001]. Hence, the crucial consideration for the ap- 
pearance of non-Fickian dispersion in a specified scale of 
a heterogeneous medium are the physical criteria for the 
power law (|3^) and its (time) range of applicability. Non- 
Fickian transport that displays these characteristics has 
been documented in several analyses of numerical sim- 
ulations, and laboratory and field data [Berkowitz and 
Scher, 1998; Hatano and Hatano, 1998; Berkowitz et al, 
2000; Kosakowski et al, 2001]. 

The large time regime of ip(s,t) corresponds to the 
small u regime for its Laplace transform and the expan- 
sion in u (for ([30j)) is quite different from (|2^) [Shlesinger, 
1974], i.e., 



xp(8,u) =p'i(s) -p' 2 (s)u p + 



(35) 



for u — > for < p < 1. Insertin g (]35 |) into (20), paral- 
lel to the development following (|26|), yields a transport 
equation from ( p^ ) which remains non-local in time and 
is not the ADE. Our development [Berkowitz and Scher, 
19951 of non-Fickian transport has been based directly 
on (|15|). In other words, solutions for the full evolution 
of a tracer plume, as well as for breakthrough curves (i.e., 
spatial and temporal distributions of tracer) can be de- 
rived directly from ( |l5| ) [e.g., Scher and Montroll, 1975; 
Berkowitz and Scher, 1997, 1998]. A (fractional) pde 
form of the transport equation, derived fro m (p2|) and 
holding only for the power law dependence (30), i.e., a 
special case of CTRW, is exhibited in section 3.2. We 
observe also that the u — > expansion of ip(s,u) for 
1 < P < 2 is similar to (|26|), but with the u 2 term re- 
placed by one proportional to u@ . In this case the correc- 
tion to the u- independent term pi(s) /fused in (p8|), j29| ) 
is proportional to u 13 ' 1 and can be significant (especially 
for/3«l). 

Finally, we note that the general CTRW formalism 
(i.e., not restricted to (|30|)) can be used to model a 
large number of physical processes. For example, ip( s ^) 
has been defined for multiple trapping [e.g., Scher et al, 
1991; Hatano and Hatano, 1998] and as such can be used 
for multiple-rate models [Haggerty and Gorelick, 1995] 
and to quantify dispersion in stratified formations [Math- 
eron and de Marsily, 1980]. Zumofen et al. [1991] have 
used the CTRW explicitly to model the latter. 



3. Fractional Differential Equations 

There is growing interest in the development and ap- 
plication of fractional differential formulations of trans- 
port equations. In particular, fractional differential equa- 



tions of the diffusion, diffusion-advection, and Fokker- 
Planck type have been considered in stochastic modeling 
in physics [e.g., Hilfer, 2000; Metzler and Klafter, 2000]. 
Here we consider fractional derivative equations (FDE) 
for transport and show how they are special cases of the 
CTRW equations developed in the previous section. We 
emphasize that FDE are not different models from the 
CTRW; rather, they are seen to emerge as asymptotic 
limit cases of the CTRW theory. 

A word of caution: referring to a transport equation 
as "fractional" can be with respect to the occurrence 
of fractional order differentiation in time or space, or 
both. Moreover, a number of definitions for fractional 
operators exist. Here, we concentrate on two possibil- 
ities: the Riemann-Liouville fractional time derivative 
o-Df (for which we will employ the more suggestive nota- 
tion d^/dt 13 ), and the Riesz spatial derivative V M [Old- 
ham and Spanier, 1974; Samko et al, 1993]. 

The development of FDE in both the time and space 
variables necessitates a more general starting equation 
than (p2|), which depends on the validity of the expan- 
sion of Pjs, t) similar to (||) . We return to the general 
solution (|2l]). In what follows, in order to obtain FDE's, 
we need the product form p(s)ifj(t) for the ijj{ s it) proba- 
bility density function, which assumes that the transition 
length and time are statistically independent quantities. 
Furthermore we need the asymptotic form ([30j) of ip(t) 
and/or p(s) (cf. below). The indicated power-law decay 
for < P < 1 causes the divergence of t, the mean tran- 
sition time (cf. section 2.4). Corresponding to ( ^p| ) the 
Laplace transform of ifj(t) is 



■ip{u) 



1 



{uc t ) 







(36) 



which is (35) summed over s, where ct is a dimen- 
sional constant determined by the physical model. Along 
the same line we consider the power-law form p(s) ~ 
c s M /|s| 1+M , < jU < 2 for the transition length, where 
c s , is analogous to c t , a dimensional constant. Similar to 
ip(t), the first and second or second (spatial) moment (s) 
of p(s) are infinite for, respectively, < /i < 1 and 
1 < fi < 2. The border case for fi = 2 is the Gaussian 
lawp(s) ~ (47rc s 2 ) -1 exp(— s 2 /(4c s )). For any symmetric 
Levy stable law p(s), the asymptotic form of the Fourier 
transform of p(s) is given by 



p(k) ~l-c^|k|" 



< fi < 2. 



(37) 



3.1 Time-FDE 



We concentrate on the case < p < 1 and /x = 2, 
for which the spatial moments are finite, but the tem- 
poral moments are infinite. We consider first the case 
with no spatial bias, £(t) = (i.e., no advective trans- 
port). Insertion of ( pq ) and the low wavenumber expres- 
sion p(k) ~ 1 — c s 2 k~Mnto ( pl| ) leads to 

1 



(38) 
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(dropping the cross term (itct)^c s 2 k 2 ) where the anoma- 
lous diffusion constant is defined as Kg = c s 2 /c t ^. The 
FDE is determined by multiplying (|3^) by the denomi- 
nator of the right side and rearranging to yield 



uV(k,u) 



1 



-^k 2 u(«- 3 P(k,u)), 



(39) 



where the dimension of the generalized diffusion constant 
is [Kg] = cm 2 sec - ^. While the two terms on the left 
correspond to dP(s,t)/dt in (s,t) space, with the ini- 
tial condition P(s, 0) = <5(s) (on both sides of ( |39| ) the 
property £{dF(t)/dt} = uF(u) - P(0) is utilized), the 
factor on the right poses the problem of finding the 
corresponding Laplace inversion. One of the definitive 
responses goes back to Riemann and Liouville who ex- 
tended the Cauchy multiple integral, in order to define 
the fractional integral, 



a-? 



P(M) 



dt' 



which possesses the important property 



P(s,t) 



P(s,u). 



(40) 



(41) 



The definition ( |40| ) explicitly includes the initial value at 
time t = 0. Note that for a negative index, /dt~@ , 
the Riemann-Liouville operator denotes fractional inte- 
gration whereas for a positive index, d 13 /dt@, we have 
fractional differentiation. In our case fractional differen- 
tiation is established as the succession of fractional inte- 
gration and standard differentiation: 



d l -p 
dt^ 



P(s,t) 



d d-P 
didt-P 



P(M). 



(42) 



With these definitions, we can now invert (39), and 
obtain the fractional diffusion equation 



r)P f) 1 -! 3 

I = ^V>P(s,<). 



(43) 



In the limit /? — * 1 ( f43[ ) reduces to the standard Brownian 
version. 

The generalization to a fractional ADE for anomalous 
transport (0 < (3 < 1), which includes a spatial bias (ad- 
vective transport), follows the same procedure as above 
[Compte, 1997; Compte et al., 1997; Compte and Caceres, 
1998; Metzler et al, 1998; Metzler and Compte, 2000], 

J^(M) = J^(-v^.V + ^V 2 )P(s,i) (44) 



where vg is the "generalized drift velocity" . Note that 
( p3| ) and (Q) involve fractional differentiation in time 
on the spatial derivative terms of the equations. These 
equations can be rewritten so they do not involve mixed 
derivatives, if desired [Metzler and Klafter, 2000]. We 
stress that the form of (E3) and (Eh relies on using (36), 



and that ([14|) is valid only for < < 1; it is modified 
significantly for 1 < f3 < 2. We have thus shown that 
the probability density P(s, f) described by the time- 
fractional ADE ([l4|) , is equivalent to the large-time limit 
of the CTRW with a bias, with the asymptotic form of 
ip(t) given by ( ^p|) (or i>(u) given by @)). For a spe- 
cific class of ip(t) (which also fulfills the asymptotic form 
(^)), the equivalence between CTRW and FDE can be 
shown over the entire range of t [Hilfer and Anton, 1995]. 



3.2. Space-FDE: Levy Flights 

We now consider the opposite case of a transition time 
distribution with an existing first moment, (3 > 1, ip(u) ~ 
1 — uct, and a transition length distribution p(s) with 
a diverging second moment (0 < fi < 2) (J-{p(s)} in 
([57])). This case can be shown to be a Markovian process 
(in contrast to the semi-Markovian process discussed in 
section 2.3) called a Levy flight. 

To avoid confusion, we stress that a Levy flight refers 
to a random movement in space, where the length of the 
transitions is considered at discrete steps, but time is not 
involved. Levy walks, on the other hand, attach a time 
"penalty" , by assigning a velocity to each transition in 
space. In the simplest case, this velocity is constant; re- 
laxation of this condition leads back to the more general 
CTRW formulation of section 2.3 [Klafter et al, 1987; 
Shlesinger et al, 1993]. In any case, Levy walks can- 
not be described in terms of simple fractional transport 
equations [Metzler, 2000]. 

A Levy flight is characterized through the Fourier- 
Laplace transform [Bouchaud and Georges, 1990; 
Compte, 1996; Metzler and Klafter, 2000] 



V(k,u) 



1 



K» k ^ 



(45) 



from which, upon Fourier and Laplace inversion, the FDE 
[Compte, 1996] 



^P( S ,()=^VP( S ,t) 



(46) 



is inferred. The Riesz operator is defined through 
[Samko et al, 1993] 



p{v^p(s,<)} = -|krp(M). 



(47) 



Note that wc use the definition K^ = c s ^ jc t for the diffu- 
sion constant. From J45|), one recovers the characteristic 
function 



P(k,i) =exp(-i^|k|' 1 ), 



(48) 



which is the characteristic function of a centered and 
symmetric Levy distribution with the asymptotic power- 
law behavior [Levy, 1925, 1954; Gnedenko and Kol- 
mogorov, 1954] 



P(s,t) - |s| 



(49) 
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Levy distributions are used to generate Levy nights 
[Bouchaud and Georges, 1990]. Accordingly, the second 
moment of a Levy flight diverges: 



«t) 2 ) 



(50) 



Observe that Levy flights are characterized by a transi- 
tion time distribution tp(t) with a finite first moment; 
they are thus fundamentally different from those pro- 
cesses underlying the time-fractional dispersion equation 
(|44|). As can be seen both descriptions are included in 
the CTRW framework. 

Including a bias into the transition distribution, one 
obtains for an asymptotic form of p(s) the Levy flight 
fractional ADE [Metzler et al, 1998] 



v • VP(s,<) = A^V At P(s,i) (51) 



which exhibits Galilei symmetry, i.e., (|5l]) is solved by 
the Levy stable solution (|49|), to be taken at the point 
s — vt. This means that the symmetric Levy stable 
plume is entirely shifted along the velocity vector v, a 
situation which strongly contrasts the growing skewness 
in the CTRW case for long-tailed transition times. Of 
course, this solution features the same divergence (^) of 
the second moment of the plume distribution. The first 
moment of (j5^) exists for all < /i < 2 and follows the 
usual Galilei symmetry expression 



(s(t)> = vt. 



3.3. Applications 



(52) 



As discussed above, although both time and space FDE 
forms are special cases of the CTRW, and both represent 
generalizations of the Fickian-based ADE, there are clear 
and critical distinctions between the transport equations 
that result from these two formulations. Here, we assess 
the Levy flight description and argue that its characteris- 
tics strongly limit its applicability to describing transport 
in geological formations. 

We consider the underlying physical picture of the 
Levy flight, as applied to tracer migration in geologi- 
cal formations: a necessary condition for the Levy flight 
description is that the domain clearly contain "streaks" 
of high and low permeability, arranged so as to lead to 
particle transitions of high and low velocity. In other 
words, the physical picture of a Levy flight requires an 
encounter with a wide range of lengths of permeability 
streaks to obtain a non-Fickian distribution of particle 
transitions. And yet, such non-Fickian distributions arise 
even without the presence of such a permeability distri- 
bution, as clearly demonstrated by, e.g., Silliman and 
Simpson [1987]. 

In addition, we observe that in mathematical terms, 
the first and second moments are often used to charac- 
terize plume migration. These quantities describe the 



spatio-temporal distribution of the tracer particles; the 
particles carry a finite mass, and therefore have a fi- 
nite velocity. As noted above, the Levy flight descrip- 
tion leads to a diverging second moment of the migrat- 
ing plume. Given that the macrodispersion parameter 
is typically defined in terms of the second moment, this 
divergence property cannot be ignored. Moreover, we 
observe that through scaling arguments [Jespersen et al., 
1999], transport only undergoes a "superdiffusive" (faster 
than linear) process; in the Levy flight description, sub- 
diffusive transport can never occur. 

With respect to the issue of a diverging second mo- 
ment, one might attempt to work with a finite number of 
sampled tracer particles in a finite range, during a finite 
time window; this leads to a truncated Levy distribution 
with finite moments. For truncated Levy distributions 
it is known that their scaling behaviors in time pertain 
up to relatively large times [Mantegna and Stanley, 1994, 
1995]. The difficulty is that to account for the temporal 
evolution of the particle cloud, the cutoffs would have to 
be adjusted to the actual space volume explored by the 
tracer particles, i.e., the cutoffs would themselves become 
time-dependent [Jespersen et al., 1999]. Put somewhat 
differently, the spatial-fractional formulation is based on 
an assumed fractal, scale- free nature of the transport pro- 
cess. Truncating the distributions leads, by definition, to 
a scale-dependent process which invalidates the use of 
simple fractional operators. 

In contrast to the above arguments, the formulation 
given by, e.g., (44), or, more generally, by ([!(])-( 19), does 
not suffer from these limitations or assumptions. In real- 
istic field situations, the distribution of particle velocities 
is expected to vary widely on the order of magnitude of 
typical spacing between sampling points. Of course, the 
velocity distribution is bounded by some maximum ve- 
locity. In the long time limit, corresponding to the small 
u limit that is of interest in our modeling, the mean effect 
of this finite variation of velocities can be approximated 
by a typical velocity. From this point of view, there- 
fore, anomalies in the plume and the related moments 
should arise from temporal "sticking" processes (i.e., low 
velocity particle transitions) which are taken into consid- 
eration in the CTRW picture. Depending on the range 
of j3 (recall (j3l"l)- (|34|)), both subdiffusive and superdiffu- 
sive behaviors for plume spreading can be characterized. 
Moreover, explicit spatial structure (well-defined conduc- 
tivity features) can be incorporated within the CTRW 
framework. 



4. The use of CTRW-based ensemble averages in 
non-stationary media: The relation of field scales 
and uncertainty 

We return now to consider the issue, raised in the In- 
troduction, that the interplay between ensemble averag- 
ing and spatial scales of non-stationary geological fea- 
tures strongly affects efforts to model transport. Broadly 
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speaking, there exist two approaches to modeling trans- 
port in large, field-scale formations. In the first ap- 
proach, the formation is treated as a single domain, 
with heterogeneities characterized and distributed ac- 
cording to a random field, with or without correlation 
and/or anisotropy. Generally speaking, these characteri- 
zations treat the domain as a stationary system, although 
stochastic models that incorporate a deterministic drift 
component (in the random field generator) have been 
considered [e.g., Li and McLaughlin, 1995]. In the second 
approach, a physical picture of the domain is constructed 
which includes explicitly specified (prescribed or known) 
heterogeneities, so that the resulting domains are non- 
stationary [e.g., LaBolle and Fogg, 2001; Koltermann and 
Gorelick, 1996; Eggleston and Rojstaczer, 1998; Feehley 
et al, 2000]. 

While the study of ensemble-averaged (stationary) do- 
mains has given rise to a sub-literature on stochastic 
methodologies and limiting behavior (e.g., perturbation 
techniques, macrodispersion) it has not yielded a prac- 
tical numerical scheme to deal with the large majority 
of field sites. Anderson [1997] describes in detail hetero- 
geneity and trending structures evident in natural geo- 
logical formations, and argues convincingly for the need 
to use facies modeling (coupled with geostatistical tech- 
niques) and/or depositional simulation models. These 
models can provide the underlying hydraulic conductivity 
structure and flow field of non-stationary domains, con- 
ditioned on field measurements, and be integrated with 
predictive models of transport. 

Within the framework of non-stationary domains, ex- 
plicitly characterized by structural trends, the question 
then arises as to how best to model transport (or, more 
precisely, how to deal with the unresolved heterogeneities 
(residues)). Clearly, there is a critical interplay between 
length scales associated with the trends and the residues. 
This gives rise to the associated uncertainty in both 
the measured/estimated hydraulic parameters and the 
measured/predicted concentrations. The generally ac- 
cepted explanation for non-Fickian transport is that het- 
erogeneities which cannot be ignored are present at all 
scales. Therefore, accounting for these residues is a cen- 
tral consideration for the quantification of non-Fickian 
transport. 

In efforts to combine non-stationarity with local-scale 
heterogeneity and uncertainty, several recent studies have 
attempted to use ADE-based modeling approaches in 
conjunction with facies modeling [e.g., Eggleston and Ro- 
jstaczer, 1998; Feehley et al., 2000]. However, these stud- 
ies, which incorporated even highly discretized systems 
(e.g., with block sizes of the order of 10 m 3 in large 
aquifers), demonstrated an inability to adequately cap- 
ture the migration patterns; these results suggests that 
unresolved heterogeneities also exist at these relatively 
small scales. We note that non-Fickian transport has 
been observed even in small-scale, relatively homoge- 
neous, laboratory-scale models [Berkowitz et al., 2000]. 
Other related issues that have been considered recently 



focus on the relative importance of diffusion and local- 
scale dispersion and on how to separate diffusive mass 
transfer processes from slow particle velocities [e.g., Har- 
vey and Gorelick, 2000; LaBolle and Fogg, 2001]. These 
questions may be considered to be somewhat moot, espe- 
cially given that "dispersion" is an artifact of averaging in 
mathematical formulations, while a definitive separation 
between diffusion and very low velocity may be unneces- 
sary. 

At all of these smaller scales, i.e., within individual fa- 
cies or depositional structures, the CTRW-based trans- 
port equations are highly effective. We therefore suggest 
that the CTRW-based approach should be used together 
with these facies and depositional models. As is usually 
done, a numerical model can be constructed which ac- 
counts explicitly for the heterogeneity structure of a for- 
mation, and the usual methods to solve for the flow field 
can be implemented. A CTRW-based transport equation 
can then be applied, rather than the ADE, over the entire 
domain. We observe that while the ADE (and the usual 
definition of "dispersion" ) is simpler to apply than the 
CTRW-based equation, the preceding discussion (both 
in this section and the previous ones) demonstrate that 
it cannot and should not generally be applied in realistic 
field situations. 

In this context, we shall consider the use of a hybrid 
model: known conductivity structures are accounted for 
explicitly, and within each block (pixel or voxel) of a 
numerical model we use the CTRW to account for the 
residues. Precluding the use of ifi(s, t) with (spatial) Levy 
forms, because the trends are included explicitly in the 
numerical model, we can start with (|2^) as a basis for 
our numerical treatment. The methods developed with 
the use of the ADE, can be carried out with the Laplace 
transforms of (E3)-(J25f), 



uP(s,u) - P (s) 



-v^(u) • VP(s,u) 
+^{u) : VVP(s,u) (53) 



uS s -0(s, u)s 
1 - ijj(u) 



uT, s ip(s, u)iss 
1 - i>(u) 



(54) 



(55) 



where Po(s) is the initial condition. 

The transport equation ( |53| ) is very similar to the 
Laplace transform of the ADE, but with the important 
exception that v^, and <f>ip are it-dependent. A spatial 
grid can be employed to numerically solve (|53"|), exactly 
as can be done with the ADE applied to a non-stationary 
system. At each grid point, the velocity value determined 
from the solution to the steady flow problem is used in 
(HMH)' a l° n g with the corresponding estimate of (3, to 
change the parameters of "0(s,it) and ip(u). 
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In this methodology the interpretation of t/>(s, u) 
changes somewhat. Instead of single transitions, we con- 
sider ip(s,u) as playing the role of accounting for the 
transition across an entire element of the spatial grid. 
This interpretation has been justified by Margolin and 
Berkowitz [2000]. 

If we insert (recall (|35D 



i)(u) ^ 1 - cpu 13 , for < (3 < 1 



(56) 



into (|53|)-(|55|), we generate non-Fickian transport across 
each block element (with cp proportional to the velocity 
value at the grid point, divided by a characteristic length, 
all raised to the (3 power). The non-Fickian behavior is 
due to the unresolved heterogeneities below the scale of 
the spatial grid. Estimates of /3 and cp can be obtained 
for each facies from a standard tracer breakthrough test 
and subsequent comparison and fitting with analytical 
solutions (as done, e.g., in Berkowitz et al. [2000] and 
Kosakowski et al. [2001]); this procedure is exactly anal- 
ogous to the usual determination of the dispersivity pa- 
rameter a in the ADE. 

Using a more complete expression for ?/>(s, u) we can 
also evolve the dynamics of the plume at very long time 
into a Gaussian (i.e., in a time regime in which ?/)(s,t) 
possesses a finite first and second temporal moment). 
The change in ^(s, u) across the boundaries can be han- 
dled by using suitable averages similar to the ADE-based 
numerical treatments. Hence one can numerically solve 
for P(s, u) at each grid point and obtain the normalized 
concentration P(s,t) by calculating £ _1 [P(s, u)]. How- 
ever, the inversion of a Laplace transform can be chal- 
lenging, and remains a key issue for future research. 

Finally, if we include pumping wells at some of the grid 
points s p (where ip(s p ,u) = 0, because the particles enter 
the well but do not emerge), then we can obtain the ac- 
cumulated concentration directly from P(s p ,u — ► 0). In 
other words, P(s p , 0) = J dtP(s p , t), and because mass 
is conserved, each pumping well acts as a sink extracting 
a fraction of the migrating particles. 



5. Summary and Conclusions 

The application of stochastic approaches to quantifica- 
tion of transport in heterogeneous geological media rests 
inevitably on the underlying conceptual picture of disper- 
sive mechanisms. The fundamental significance of this 
picture was pointed out long ago. As noted by Bear 
[1972], in his discussion of the work of Scheidegger [1954, 
1958], "...the application of the statistical approach re- 
quires... a choice of the type of statistics to be employed, 
i.e., the probability of occurrence of events during small 
time intervals within the chosen ensemble. This may 
take the form of correlation functions between velocities 
at different points or different times, or joint-probability 
densities of the local velocity components of the parti- 
cle as functions of time and space or a probability of an 



elementary particle displacment. The chosen correlation 
function determines the type of dispersion equation de- 
rived." 

We have developed this early insight into a full, quan- 
titative theory where the joint probability density is the 
^(s, t). This joint spatial-temporal distribution allows us 
to account for the behavior of migrating particles which 
can encounter a wide range of velocity regions in het- 
erogeneity lenses of different spatial dimensions. This 
approach is in contrast to most others which have, his- 
torically, emphasized spatial formulations of transport 
equations, motivated by the clear spatial heterogeneity 
of geological formations. 

The overarching framework for our physical picture of 
transport, and the assumptions (as detailed above) on 
particle transitions, is the Master Equation. This equa- 
tion represents a general, yet highly applicable, quantifi- 
cation of transport which recognizes the broad spectrum 
of particle motions in space and time. We show, under 
a general assumption of the form of w(s, s'), that the 
Master Equation can be specialized in any single realiza- 
tion of the geological domain to a generalized form of the 
ADE. 

The ensemble average of the unrestricted Master Equa- 
tion leads to a Generalized Master Equation, which is 
exactly equivalent to the CTRW. As a limiting form, un- 
der highly restrictive conditions regarding the character 
of the transport (and therefore of the degree of structural 
heterogeneity), the conventional ADE can be recovered 
from this formulation. 

Aquifers are inherently heterogeneous over a wide 
range of scales, and Fickian transport (embodied in the 
ADE) does not generally occur on practical scales of in- 
terest. We therefore suggest that the overwhelming fo- 
cus on defining "effective" dispersion, or "macrodisper- 
sion" coefficients, in Fickian or pseudo-Fickian formula- 
tions of the transport problem, is misplaced for field-scale 
problems. The CTRW theory, which is the basis for our 
transport equation, quantifies naturally the non-Fickian 
behavior observed at laboratory and field scales, as well 
as in numerical simulations. The essential character of 
the transport can be embodied in an asymptotic form of 
the ip(s, t), specifically by an exponent (3. This exponent, 
which can be determined from the velocity distribution 
(based on solution of flow for a given conductivity field) 
or from a tracer test, parameterizes an entire class of 
non-Fickian plume evolutions, on scales larger than the 
size of the heterogeneities. Detailed discussions on the 
practical identification of -0(s,i) and parameter values is 
given in Berkowitz and Scher [2001], Kosakowski et al. 
[2000], and Berkowitz et al. [2000, 2001]. 

We have also shown how fractional derivative formu- 
lations of transport equations are special, asymptotic 
(limit) cases, (§0|) for ip(sJ), of the CTRW theory. In- 
serting this limiting form ( |35| ) into the Laplace transform 
of (|l^), one arrives at the same step necessarily encoun- 
tered at the outset of the solution of the FDE. Retention 
of the more general equation (|l^) has important advan- 
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tages for a more complete modeling of the transport pro- 
cess. The limiting forms characterized by the exponent 
(3 (which is the fractional order of the derivative in the 
FDE) apply for a certain time range only. Beyond this 
range, the rp(s,t) changes in a manner that allows the 
plume to eventually assume a Gaussian shape (defined 
by "macrodispersion" ) as is reasonable for most physical 
systems. 

Finally, we consider how best to quantify contaminant 
transport in non-stationary geological formations. We 
delineate a hybrid approach in which known structural 
properties are included explicitly, and unresolved (un- 
known) heterogeneities at smaller scales are accounted for 
within the CTRW theory. Practical application of this 
approach is achieved by replacing the usual ADE equa- 
tion that is integrated into numerical simulation codes 



by a CTRW-based transport equation. This transport 
model can be integrated with existing numerical model- 
ing techniques to determine the underlying flow field. 

We are currently focusing efforts on implementation 
of the solution technique suggested here, as well as on 
deriving analytical solutions for CTRW-based transport 
equations for forms of ip(s, t) generalized in both space 
and time. 

Appendix A 

We showed how the use of (||)-(|8|) leads to the expres- 
sion ( |l0| ) for the first term of the right side of ||). We 
outline the derivation here for the second and third terms 
of the right side of (^), using these same equations. We 
have for the second term 



^w(s,s')(s'-s).VC(s,() « ^(F(|s'-s|;s) + (s'-s)- VF) [A + ±(s' - s)Vtt)] (s'-s)-VC(s,t) 

« ]T AF(|s' - s|; s)i(s' - s)^V - s) • VC(s, *) 

s' 

+ ^(s' - s) • VFA(s' - s) • VC(s, t) 



D(s)^VC(s, t) + 2VD(s)VC(s, t) 
A 



(57) 



and for the third term, 

^iu(s,s')|(s'-s)(s'-s) : VVC(s,t) «E S ' F (I S ' -s|;s)A±(s' -s)(s' - s) : VVC(s,t) = D(s) : VVC(s,t) (58) 

s' 

I 



We add the results of (Al), (A2) and (jig) to obtain 
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